A fast and accurate algorithm for solving 
Bernstein— Vandermonde linear systems 



Ana Marco 1 , Jose- Javier Martinez 2 

Departamento de Matemdticas, Universidad de Alcald, 
Campus Universitario, 28871-Alcald de Henares (Madrid), Spain 



Abstract 

A fast and accurate algorithm for solving a Bernstein- Vandermonde linear sys- 
tem is presented. The algorithm is derived by using results related to the bidiagonal 
decomposition of the inverse of a totally positive matrix by means of Neville elimi- 
nation. The use of explicit expressions for the determinants involved in the process 
serves to make the algorithm both fast and accurate. 
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1 Introduction 



The Bernstein basis for the space of algebraic polynomials of degree less than 
or equal to n is a widely used basis in Computer Aided Geometric Design due 
to the good properties that it possesses (see, for instance, [5, 9, 10, 11, 17]). 
However, the explicit conversion between the Bernstein and the power basis is 
exponentially ill-conditioned as the polynomial degree increases [10]. For this 
reason, it is very important that when designing algorithms for performing 
numerical computations with polynomials expressed in Bernstein form, all 
the intermediate operations are developed using this form only [2]. A paper 
which presents various basic operations for polynomials in Bernstein form is 
[11]. In [2] an algorithm for computing the greatest common divisor of two 
polynomials in Bernstein form avoiding explicit basis conversion is given. 
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Our aim in this paper is to develop a fast and accurate algorithm for solving 
a linear system whose coefficient matrix is a Bernstein- Vandermonde matrix, 
that is to say, a Vandermonde-like matrix for the Bernstein polynomials. Tak- 
ing into account that this matrix is the coefficient matrix of the linear system 
associated with a Lagrange interpolation problem in the Bernstein basis, this 
result will allow us to perform a basic polynomial procedure, the Lagrange in- 
terpolation [7], avoiding transformations between Bernstein and power basis. 
In this way, general algorithms for polynomials in Bernstein form involving an 
interpolation stage could be designed using this form only. 

Our algorithm will be based on the bidiagonal factorization of the inverse 
of the Bernstein- Vandermonde matrix. Factorizations in terms of bidiagonal 
matrices are very useful when working with Vandermonde [3, 14, 16], Cauchy 
[4], Cauchy- Vandermonde [20, 21] and generalized Vandermonde matrices [8]. 

Let us observe here that, of course, Bernstein- Vandermonde linear systems 
can be solved by using standard algorithms such as Gaussian or Neville elim- 
ination. However they are not fast and the solution provided by them will 
generally be less accurate since Bernstein- Vandermonde matrices are ill con- 
ditioned (see Table 1, where n is the degree of the Bernstein polynomials 
involved in the definition of the matrix, whose order is n + 1) and these algo- 
rithms (which do not take into account the structure of the matrix) can suffer 
from subtractive cancellation. 

Table 1 

Condition number of a Bernstein- Vandermonde matrix with Xj = — ^ . 



n 


10 


15 


20 


25 


30 


35 


40 


K>oo (A n ) 


2.1e + 04 


2.6e + 06 


3.5e + 08 


4.7e + 10 


6.6e + 12 


9.0e+ 14 


1.3e + 17 



On the other hand, it must be observed that the high condition numbers are 
due to the high norm of the inverse matrix, since as it will be easily seen after 
the definitions of Section 3 the oo-norm of a Bernstein- Vandermonde matrix is 
always 1. Therefore, taking into account that if Ax = b and A(x+Ax) = b+Ab 
then Ax = A _1 Ab, the fact that 1 1 ^4 1 j | is high implies that the effect of 
perturbations in the vector b is likely to be important. A lot information 
about the related concepts of perturbation theory and numerical stability of 
algorithms in the context of solving linear systems can be found in Chapter 7 
of [16]. 

The rest of the paper is organized as follows. Neville elimination and total pos- 
itivity are considered in Section 2. In Section 3 the bidiagonal factorization of 
the inverse of a Bernstein- Vandermonde matrix is presented. In Section 4 the 
algorithm for solving a linear system whose coefficient matrix is Bernstein- 
Vandermonde, and the computation of its complexity are given. Finally, Sec- 
tion 5 is devoted to illustrate the accuracy of the algorithm by means of some 
numerical experiments. 
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2 Basic results on Neville elimination and total positivity 



In this section we will briefly recall some basic results on Neville elimination 
and total positivity which we will apply in Section 3. Our notation follows the 
notation used in [12] and [13]. Given k, n 6 N (1 < K n), Qk,n will denote 
the set of all increasing sequences of k natural numbers less than or equal to 
n. 

Let A be a real square matrix of order n. For k < n, m < n, and for any 
ot £ Qk,n and j3 G Q m ,n, we will denote by A[a|/3] the submatrix k x m of A 
containing the rows numbered by a and the columns numbered by (3. 

The fundamental tool for obtaining the results presented in this paper is the 
Neville elimination [12, 13], a procedure that makes zeros in a matrix adding 
to a given row an appropriate multiple of the previous one. For a nonsingular 
matrix A = (ai,j)i<i,j< n , it consists on n — 1 steps resulting in a sequence of 
matrices A := Ai — > A 2 — > . . . — > A n , where A t = {o%j)i<i,j< n has zeros below 
its main diagonal in the t — 1 first columns. The matrix A t+1 is obtained from 
A t (t — 1, . . . , n) by using the following formula: 



(t+i) 



«S - (4VaU*K-U if ^ > t + 1 and j > t + 1 . (2.1) 
otherwise 



a^ j if i < t 



In this process the element 

Pij := alj 1 < j < n; j <i<n 

is called pivoi of the Neville elimination of A. The process would break 
down if any of the pivots p^j (j < i < n) is zero. In that case we can move 
the corresponding rows to the bottom and proceed with the new matrix, as 
described in [12]. The Neville elimination can be done without row exchanges 
if all the pivots are nonzero. The pivots p^ are called diagonal pivots. If all 
the pivots Pij are nonzero, then p^ = a it i Vi and, by Lemma 2.6 of [12] 

det A\i — j + 1, . . . , ill, . . . , il 

^ 51R^CT 1<J - ! -"' (2 ' 2) 



The element 

m id = 1<] <n; j <i<n (2.3) 

Pi-l,j 

is called multiplier of the Neville elimination of A. The matrix [/ := A n is 
upper triangular and has the diagonal pivots in its main diagonal. 
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The complete Neville elimination of a matrix A consists on performing the 
Neville elimination of A for obtaining U and then continue with the Neville 
elimination of U T . The pivot (respectively, multiplier) of the complete 

Neville elimination of A is the pivot (respectively, multiplier) (j, i) of the 
Neville elimination of U T , if j > i. When no row exchanges are needed in the 
Neville elimination of A and U T , we say that the complete Neville elimination 
of A can be done without row and column exchanges, and in this case the 
multipliers of the complete Neville elimination of A are the multipliers of the 
Neville elimination of A if i > j and the multipliers of the Neville elimination 
of A T if j > i. 

A matrix is called totally positive (respectively, strictly totally positive) if all 
its minors are nonnegative (respectively, positive). The Neville elimination 
characterizes the strictly totally positive matrices as follows [12]: 

Theorem 2.1. A matrix is strictly totally positive if and only if its complete 
Neville elimination can be performed without row and column exchanges, 
the multipliers of the Neville elimination of A and A T are positive, and the 
diagonal pivots of the Neville elimination of A are positive. 

It is well known [5] that the Bernstein- Vandermonde matrix is a strictly totally 
positive matrix when the interpolation points satisfy < x\ < x 2 < . . . < 
x n+ i < 1, and this fact has inspired our search for a fast algorithm, but this 
result will also be shown to be a consequence of our Theorem 3.3. 



3 Bidiagonal factorization 

The Bernstein basis of the space H n (x) of polynomials of degree less than or 
equal to n on the interval [0, 1] is: 




We will call the following Vandermonde-like matrix for the Bernstein basis B n , 

A= ©a-^r (?> 2 (i - ••• {i)x$ 
v (o)(i-^ + i) n GK+ia-^+i)"- 1 ••• (:;)< + J 

a Bernstein- Vandermonde matrix. Let us observe that these Vandermonde- 
like matrices are not exactly the Vandermonde-like matrices considered in 
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[16], since all the polynomials in the Bernstein basis B n have the same degree 
n and they do not satisfy a three-term recurrence relation. From now on, we 
will assume < x 1 < x 2 < . . . < x n+1 < 1. 

This matrix A is the coefficient matrix of the linear system associated with the 
following Lagrange interpolation problem in the Bernstein basis B n : given the 
interpolation nodes {x{ : i — 1, . . . ,n + l} with < X\ < x 2 < ■ ■ ■ < x n+ \ < 1 
and the interpolation data {bi : i — 1, . . . , n + 1} find the polynomial 



p(x) = ^a k ^j(l-x) 



n-k^k 



such that p(xi) =b; t for i — 1, . . . , n + 1. A good introduction to the interpo- 
lation theory can be seen in [7]. 



Proposition 3.1. 



det A 



n 
,0, 



,1, 



n' 



n 



Xi). 



l<i<j<n+l 



Proof. It is easy to see that the matrix of change of basis from the Bernstein 
basis B n to the power basis {l,x,x 2 , . . . ,x n } is a lower triangular matrix of 
order n+1 whose diagonal elements are (o) > (") > • • • > (") • 



From this fact, it is obtained that 

det A 



n \ In 
0/U, 



n 1 



det V, 



where V is the Vandermonde matrix 



V 



1 X\ x\ 
1 X2 x\ 



n \ 



\l x n+ i x\ +l ■ ■ ■ x™ +1 j 
Using the well-known formula for the determinant of a Vandermonde matrix 

det V — | j (xj — x^ 

l<i<j<n+l 



the proof is concluded. □ 



The next corollary follows directly from Proposition 3.1., and will be useful to 
make the derivation of the algorithm easier. 
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Corollary 3.2. 



det 



(1- Xl ) n x^l-x^ 1 - 1 
(l-x 2 ) n x 2 (l-x 2 ) n - 1 



X 



X 



\{l - X n+ i) n X n+ i(\ - X n+ i) 



n-l 



X 



n 



(Xi 



Xi) 



l<i<j<n+l 



n+1/ 



The following result will be the key to construct our algorithm. 



Theorem 3.3. Let A = (oj,j)i<jj<n+i be a Bernstein- Vandermonde matrix 
whose nodes satisfy < x 1 < x 2 < . . . < x n < x n+ i < 1. Then A^ 1 admits a 
factorization in the form 



A 1 — G\G 2 - • • G n D 1 F n F n -i ■ ■ ■ F±, 



(3-1) 



where Gi are upper triangular bidiagonal matrices, Fj are lower triangular 
bidiagonal matrices (i — 1, . . . , n), and D is a diagonal matrix. 

Proof. The matrix A is strictly totally positive (see [5]) and therefore, by 
Theorem 2.1, the complete Neville elimination of A can be performed without 
row and column exchanges providing the following factorization of A~ Y (see 
[12] and [13]): 

A 1 = G\G 2 ---G n D 1 F n F n _i • • • Fi, 
where Fi (1 < % < n) are bidiagonal matrices of the form 



1 

1 



1 

-rrii+i,i 1 

-mi+2,i 1 



-m n+lji 1 



(3.2) 



6 



Gj (1 < % < n) are bidiagonal matrices of the form 



Gj = 



1 

1 



1 

-rn i+1>i 1 

-rn i+2 ,i 1 



V 



-m n+lji 1 



(3.3) 



and D is the diagonal matrix whose ith (1 < % < n + 1) diagonal entry is the 
diagonal pivot p ifi = afj of the Neville elimination of A: 



D = diag{pi 5l ,p 2)2 . . . ,p„+i,„+i}. 



(3.4) 



Taking into account that the minors of A with j initial consecutive columns 
and j consecutive rows starting with row % are 



det A[i, + j - 



;)(")- (A) 

i<fc</<i+j-l 



a result that follows from the properties of the determinants and Corollary 3.2, 
and that m^j are the multipliers of the Neville elimination of A, we obtain 
that 

_ (i - x t r^ + \i - Xi-j) nUfa - x t . k ) 



m i:j = 



Pi 



Pi-u ~ (l-^-o^-^nLsfa-i-^-fc) 

where j — 1, . . . , n and i — j + 1, . . . , n + 1. 



(3.5) 



As for the minors of A T with j initial consecutive columns and j consecutive 
rows starting with row i, they are: 



det A T [i,... : i + j-l\l : ...J} 



n \ ln\ I n 
K i-l)[i)'"[i + j-2J- L 

(l- Xl )»-i-i+ 2 (l-x 2 ) n - i -* +2 ---(l-x j ) n - i -* +2 1] fa -a*). 

i<fc<Kj 



i— 1 i— 1 



i-1 



This expression also follows from the properties of the determinants and Corol- 
lary 3.2. Since the entries m^j are the multipliers of the Neville elimination of 
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using the previous expression for the minors of A T with initial consecutive 
columns and consecutive rows, it is obtained that 



On — i + 2) • Xj 

mi 'i = T — V j = 1, . . . , n; % = j + 1, . . . , n + 1. (3.6) 

Finally, the ith diagonal element of D 



PM = - — ^-i z; -i t = l,...,n+l, (3.7) 



is obtained by using the expression for the minors of A with initial consecutive 
columns and initial consecutive rows. □ 

Moreover, by using the same arguments of [20], it can be seen that this factor- 
ization is unique among factorizations of this type, that is to say, factorizations 
in which the matrices involved have the properties shown by formulae (3.2), 
(3.3) and (3.4). 



Let us observe that the formulae obtained in the proof of Theorem 3.3 for 
the minors of A with j initial consecutive columns and j consecutive rows, 
and for the minors of A T with j initial consecutive columns and j consecutive 
rows show that they are not zero and so, the complete Neville elimination of 
A can be performed without row and column exchanges. Looking at equations 
(3.5)-(3.7) is easily seen that mij, rhij and are positive. Therefore, taking 
into account Theorem 2.1, this confirms that the matrix A is strictly totally 
positive. 



4 The algorithm 



In this section we will present a fast algorithm for solving a linear system whose 
coefficient matrix is a Bernstein- Vandermonde matrix. In order to solve the 
linear system Ax = b, where A is the (n + 1) x (n + 1) Bernstein- Vandermonde 
matrix introduced in Section 3, we use Theorem 3.3 for obtaining 

x = A~ x b = G l G 2 ■ ■ ■ G n D~ x F n F n _i ■ ■ ■ F 1 b. 



Since Fj and Gj (i — 1, . . . , n+1) are bidiagonal matrices and D' 1 is a diagonal 
matrix, it is clear that the computational complexity of computing the whole 
product from right to left is 0(n 2 ). It remains to see that the construction 
of the matrices F iy G; L and D^ 1 can be carried out with a computational 
complexity of 0(n 2 ). 
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Let us start with the entries m^j given by equation (3.5). We will use the 
following expressions 



m iA = (1 - Xi-i) ■ M iA 

(l-x i ) n -^ 1 Ui=i(xi-x i . k ) 



(i - ^-i)"-^ +2 ni =2 (^-i - x^u) 

= (l-Xi JiXi-Xi-j) 

'■' +J (l-xO^-i-Xi-j-i) J 
"lij+i = (1 - Zi-j-i) • M iJ+ i, 



where i — 2, . . . , n + 1 and j = 1, . . . , i — 2, in their construction: 
for i = 2 : n + 1 

m iA = (1 - Xi-i) • M 
f or j = 1 : % — 2 

(l-Xi)(Xj_l-Xi_j_l) 

m iJ+ i = (1 - Xi-j-i) ■ M 
end 
end 

Now we compute the entries m it j given by equation (3.6): 
for j — 1 : n 



for i— j + n + 1 

end 
end 

As for the diagonal elements p^j of D given by equation (3.7), they are con- 
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structed using the equalities 



Pi,i 



1 




nrii(i-^) 



n — i + 1 



Pi+l,i+l 



Qi+l,i+l 




k<i+l 



where i = 



1, . . . , n, in the following way: 



q 



= l 



pi,i = (i 




n 



for i = 1 : n 



aux = 1 
for k = 1 : i 

aux = (x i+ i — Xk) ■ aux 
end 

= g • (1 - Xi+i)™-* • crax 



Looking at this algorithm is enough to conclude that the computational com- 
plexity of the construction of the matrices D, Fi and Gi (i = 1, . . . ,n + 1) 
is 0(n 2 ), and therefore, the computational complexity of solving the whole 
linear system is also 0{n 2 ). 

A similar algorithm with computational complexity 0{n 2 ) can be developed 
for computing the bidiagonal factorization of the inverse of A T , that is, A~ T . In 
consequence, an algorithm with computational complexity 0(n 2 ) is obtained 
for solving the dual linear system A T x = b, where A is an (n + 1) x {n + 1) 
Bernstein- Vandermonde matrix, by using: 



x = A~ T b = FfFZ ■ ■ ■ F^GlGU ■ ■ -G\. 



end 
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5 Numerical experiments and final remarks 

Finally we present some numerical experiments which illustrate the good prop- 
erties of our algorithm. We compute the exact solution x e of each one of 
the Bernstein- Vandermonde linear systems Ax = b by using the command 
linsolve of Maple 10 and use it for comparing the accuracy of the results 
obtained in Matlab by means of: 

(1) The algorithm presented in Section 4 for computing the bidiagonal de- 
composition of A' 1 . We will call it MM. 

(2) The algorithm TNBD of Plamen Koev [18] that computes the bidiagonal 
decomposition of A" 1 without taking into account the structure of A. 

(3) The command A\b of Matlab. 

The fast product (from right to left) of the bidiagonal matrices and the vector 
b is also implemented in Matlab and is the second stage in the solution of the 
linear system in (1) and (2). 

We compute the relative error of a solution x of the linear system Ax = b by 
means of the formula: 

\\ x ~ x e\\2 

err = — n — n . 

\\ x eh 

Remark. The algorithm TNBD computes the matrix denoted as BV(A) in [19], 
which represents the bidiagonal decomposition of A. But it is a remarkable 
fact that the same matrix BT>(A) also serves to represent the bidiagonal de- 
composition of A -1 . The algorithm computes BT>(A) by performing Neville 
elimination on A, which involves true substractions, and therefore does not 
guarantee high relative accuracy. 

A detailed error analysis of Neville elimination, which shows the advantages 
of this type of elimination for the class of totally positive matrices, has been 
carried out in [1] , and related work for the case of Vandermonde linear systems 
can be seen in Chapter 22 of [16]. 

Example 5.1 Let £>io be the Bernstein basis of the space of polynomials with 
degree less than or equal to 10 on [0, 1] and A be the Bernstein- Vandermonde 
matrix of order 11 generated by the nodes x { = for i = 1, . . . , 11. The 
condition number of A is k 2 (A) = 1.8e + 04. Let us consider 

b\ = (1, 0, 2, -1, 3, 1, -2, 0, 0, 3, 5) T 
b T 2 = (1,-2, 1,-1,3, -1,2, -1,4, -1,1) T 

two vectors of data. 
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The relative errors obtained when using the approaches (1), (2) and (3) for 
solving the systems Ax = bi (i — 1,2) are reported in Table 2. 



h 


MM 


TNBD 


A\h 


h 


1.3e-15 


7.8e-14 


5.4e-14 


bi 


8.6e-16 


8.2e-14 


1.0e-14 



Table 2 

Example 5.1: A is a Bernstein- Vandermonde matrix of order 11. 

The following example will show how the accuracy of our approach is main- 
tained when the order, and therefore the condition number (see Table 1), of 
the Bernstein- Vandermonde matrix increases, while the accuracy of the other 
two approaches which do not exploit the structure of the matrix A goes down. 

Example 5.2 Let £>i 5 be the Bernstein basis of the space of polynomials with 
degree less than or equal to 15 on [0, 1] and A be the Bernstein- Vandermonde 
matrix of order 16 generated by the nodes Xi = -fj for % = 1, . . . , 16. The 
condition number of A is k 2 (A) = 2.3e + 06. Let 

b\ = (2, 1, 2, 3, -1, 0, 1, -2, 4, 1, 1, -3, 0, -1, -1, 2) T 

b\ = (1, -2, 1, -1, 3, -1, 2, -1, 4, -1, 2, -1, 1, -3, 1, -4) T 

two vectors containing the data. The relative errors of the solutions of the 
linear systems Ax = bi (i = 1, 2) obtained by means of the approaches (1), (2) 
and (3) are reported in Table 3. 



h 


MM 


TNBD 


A\h 


h 


1.0e-15 


5.9e-ll 


6.5e-12 


bi 


4.9e-16 


5.9e-ll 


6.4e-12 



Table 3 

Example 5.2: A is a Bernstein- Vandermonde matrix of order 16. 

Let us observe that in the first stage of the algorithm, which corresponds to the 
computation of the bidiagonal decomposition of A -1 , the high relative accuracy 
of the algorithm is obtained because no subtractive cancellation occurs: we are 
multiplying, dividing, or adding quantities with the same sign, or forming 1—Xi 
and Xi ± xj where Xi and initial data. 

As for the second stage of the algorithm corresponding to the evaluation of 
the product 

GiG 2 ■ ■ ■ G n D F n F n _i ■ ■ ■ Fib, 

it must be noted that when the vector b has alternating sign pattern ( sign(bi) = 
(±1)*), this stage is substraction-free, and so the product can be computed 
with high relative accuracy. This is a consequence to the checkerboard sign 
pattern of A -1 , which derives from the fact that A is a totally positive matrix. 
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This important property was already observed in an analogous situation in 
the paper [15], devoted to the error analysis of the Bjorck-Pereyra algorithm 
for Vandermonde systems. There the use of the increasing ordering for the 
interpolation points is recommended, an ordering which in the case < X\ < 
x 2 < ■ ■ ■ < x n makes the Vandermonde matrix a totally positive one. 

Our next two examples will serve to illustrate the concept of effective well- 
conditioning introduced by Chan and Foulser in [6]. This concept has also 
been studied in the context of totally positive Cauchy systems in [4], where 
it is seen that the BKO algorithm exploits the effective well-conditioning to 
produce higher accuracy for special right-hand sides. 

Let A e R nxn with singular value decomposition A = UY>V T , where £ = 
diag(o-j), <7i > a 2 > • • ■ > a n and U = [u\U 2 • ■ • u n ]. Let = U^U]: ', with U k = 
[u n+ i^k ■ ■ - Un], the projection operator onto the linear span of the smallest k 
left singular vectors of A. The Chan-Foulser number [4, 6] for the linear system 
Ax = f is defined by 

l A t\ ■ a n-k+l H/H2 

7(A, /) = mm — . 

k o- n \\Pkjh 

Let us observe that the computation of the Chan-Foulser number simplifies 
when / = Ui (i — 1, . . . , n). In this case, using the fact that the left singular 
vectors u±, . . . , u n are the columns of the orthogonal matrix U, we obtain 

-f(A,Ui) = — . 

Example 5.3 Let A be the same matrix as the one considered in Example 
5.2, whose condition number is k 2 (A) = 2.3e + 06. We have solved the 16 
linear systems Ax = Ui, where u±, . . . , uie are the left singular vectors of A. 
The results obtained when using the approaches (1), (2) and (3), and the 
corresponding Chan-Foulser numbers are reported in Table 4. 

Example 5.4 Let B15 be the Bernstein basis of the space of polynomials with 
degree less than or equal to 15 on [0, 1] and A be the Bernstein- Vandermonde 
matrix of order 16 generated by the interpolation nodes 

1 1 1 1 1 1 1 1 11 19 17 15 11 9 7 5 
— < — < — < — < — <-<-<-< — < — < — < — < — < — < — <-. 

18 16 14 12 10 8 6 4 20 34 30 26 18 14 10 6 

Its condition number is k 2 (A) = 3.5e + 09. We have solved the 16 linear 
systems Ax = u^, where u±, . . . , Uiq are the left singular vectors of A. 

The results obtained when using the approaches (1), (2) and (3), and the 
corresponding Chan-Foulser numbers are reported in Table 5. 
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Ui 






1 1\ JD U 




U\ 


z.oe-|-uo 


1 1 1 n 

1.16-1U 


A 7c 19 

4. / e-iz 


9 80 1 1 
z.06-11 


U2 


z. ie-j-uo 


C Od 11 

o.ue-11 


1 (Id 11 


1 7„ 11 

1. 1 e-11 


u 3 


1. 1 ctUD 


9 c„ 11 




1.0e-lU 


W4 


i.oe+uo 


/I Qo 1 1 

4.ye-n 


Q go 11 


0. 46-11 


U5 


y .oe-ruo 


4.oe-ll 


1 1 1 1 

1.16-11 


1 Ao 1 1 
1.46-11 


UQ 


D.U6-|-U0 


Q 1„ 11 

o.ie-11 


1 7„ 11 

1 . / e- 1 1 


fi 9o 1 9 
D.Z6-1Z 


Uj 


O.06-|-U0 


A Ho 1 1 
4.Ue-ll 


7 c„ 19 

/ .oe-iz 


9 1 1 1 
Z. 16-11 


Ug 


I.octUO 


1 80 1 9 


1 80 1 1 

l.oc- 1 1 


1 /Lo 1 9 
1.46-1Z 


Ug 


a c;„_i_n/i 

O.06-rU4 


1 Oo 11 

1 . ze- 1 1 


Q So 1 9 

y.oe-iz 


D.U6-1Z 


UlO 


o.4e-rU4 


1 7o 1 9 

1. / e-iz 


1 c„ 11 

1 . oe- 1 1 


Q 1 1 9 

y. i6-iz 


U\i 


l.ze-rU4 


4.ye-lo 


9 Ho 1 1 


Q 1„ 19 

0. 16-1 Z 


U\2 


o.4e-rUo 


^0 1 


1 /lo 1 1 

1.4:6-11 


1 /lo 1 1 

1.46-11 


^13 


7 Qp_i_n9 


1 zip 1 ^ 


9 4p 1 1 

Z . 4:6- 1 1 


9 4p 1 1 
Z.46-1 1 


lt!4 


1.4e+02 


8.1e-14 


6.3e-12 


2.4e-ll 


Ulb 


1.6e+01 


7.1e-15 


2.1e-ll 


1.3e-ll 




1 


5.1e-16 


5.9e-ll 


6.3e-12 



Table 4 

Example 5.3: left singular vectors Ui for the right hand side. 

The results appearing in Table 4 and Table 5 show that, as it could be 
expected, the solution vectors obtained by conventional methods based on 
Neville elimination or on Gaussian elimination suffer from a loss of digits of 
precision which is roughly proportional to the decimal logarithm of the con- 
dition number of the Bernstein- Vandermonde matrix. 

On the contrary, the results in Table 4 and Table 5 indicate that the accu- 
racy obtained by our algorithm, which exploits the structure of the matrix, 
is governed by the Chan-Foulser number. In particular, the very high relative 
precision observed in the 16th system of both examples reflects two facts: the 
Chan-Foulser number is equal to 1 and the vector w 16 has alternating sign 
pattern ( signipi) = (±1) 1 ). 
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1 1 ns 
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9 9o ns 
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u 13 




9 4p 1 9 
z . *±c - 1 z 


1 ns 


^ ns 

0. oe-uo 


Uu 


9.3e+02 


7.3e-13 


2.8e-10 


1.9e-08 


U15 


1.4e+01 


4.3e-14 


3.0e-09 


3.6e-08 


U16 


1 


7.6e-15 


l.le-09 


1.4e-08 



Table 5 

Example 5.4: left singular vectors n« for the right hand side. 
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